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Abstract 

The Vicsek model is a very popular individual based model which de- 
scribes collective behavior among animal societies. A macroscopic version 
of the Vicsek model has been derived from a large scale limit of this indi- 
vidual based model [13]. In this work, we want to numerically validate this 
Macroscopic Vicsek model (MV). To this aim, we compare the simulations 
of the macroscopic and microscopic models one with each other. The MV 
model is a non-conservative hyperbolic equation with a geometric constraint. 
Due to the lack of theory for this kind of equations, we derive several equiv- 
alents for this system leading to specific numerical schemes. The numerical 
simulations reveal that the microscopic and macroscopic models are in good 
agreement provided that we choose one of the proposed formulations based 
on a relaxation of the geometric constraint. This confirms the relevance of the 
macroscopic equation but it also calls for a better theoretical understanding 
of this type of equations. 
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1 Introduction 



This paper is devoted to the numerical study of a macroscopic version of the Vic- 
sek model which describes swarming behavior. This macroscopic model has been 
derived in [13] from the microscopic Vicsek model [29]. The goal of this work is 
to provide a numerical validation of the macroscopic model by comparing it with 
simulations of the microscopic model. 

The Vicsek model [29] is widely used to describe swarming behavior such as flock 
of birds [3], schools of fish [2,10,20,26] (in this case the model is combined with 
an attractive-repulsive force) or recently the motion of locusts [6]. In this model, 
individuals have a constant velocity and they tend to align with their neighbors. 
Despite the simplicity of the model, a lot of questions remain open about it. A 
first field of research concerns phase transitions within the model depending on the 
level of noise [7,17,25,29]. Another question arises from the long time dynamics 
of the model [11, 12, 18]: is there convergence to a stationary state of the system? 
From another perspective, since collective displacements in natural environment 
can concern up to several million individuals, it is natural to look for a macroscopic 
version of the Vicsek model. On the one hand, macroscopic models constitute 
powerful analytical tools to study the dynamics at large scales [9, 14,24]. On the 
other hand, the related numerical schemes are computationally much more efficient 
compared with particle simulations of a large number of interacting agents. In [13], 
a Macroscopic Vicsek model (MV) has been derived from a large scale limit of 
the microscopic Vicsek model. The macroscopic model is obtained from a rigorous 
perturbation theory of the original Vicsek model. Another macroscopic model is 
obtained in [4] based on more phenomenological closure assumptions. 

The MV model presents several specificities which make the model interesting. 
First, it is a non-conservative hyperbolic system and secondly it involves a geometric 
constraint. These are the consequences at the macroscopic level of two specificities 
of the microscopic model: the total momentum is not conserved by the particle 
dynamics and the speed of the particles is constant. The first property is an intrinsic 
property of self-propelled particles and the second property is a usual assumption in 
the models of collective displacements [10,16,29]. Up to our knowledge the theory 
of such systems is almost empty. Non-conservative systems have been studied in 
the literature [5,8,21] but none of them involve geometric constraints. 

In this work, since a theoretical framework for such systems is not available, we 
adopt several approaches. First, we introduce a conservative formulation for the 
ID formulation of the MV model which is equivalent to the initial one for smooth 
solutions only. With this conservative formulation, we can use standard hyper- 
bolic theory to build Riemann problem solution and shock capturing schemes [22]. 
The numerical scheme based on the conservative formulation is called conservative 
method. But since the equivalence with the original formulation is only valid for 
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smooth solutions [23] , there is no guarantee that the conservative formulation gives 
the right answer at shocks. For this reason, we introduce another formulation of 
the MV model where the constraint is treated through the relaxation limit of an 
unconstrained conservative system. This formulation leads to a natural numerical 
scheme based on a splitting between the conservative part of the equation and the 
relaxation. This scheme will be referred to as the splitting method. For compari- 
son purposes, two other numerical schemes are also used, an upwind scheme and a 
semi- conservative one (where only the mass conservation equation is treated in a 
conservative way). 

The numerical simulations of the MV model reveal that the numerical schemes 
all agree on rarefaction waves but disagree on shock waves. To determine the correct 
solution, we use the microscopic model in a regime where its solution is close to 
that of the macroscopic model. In practice, this corresponds to regimes where the 
number of particles per domain of interaction is high. The splitting method turns 
out to be in good agreement with particle simulations of the microscopic model, 
by contrast with the other schemes. In particular, for an initial condition with 
a contact discontinuity, the solution given by the conservative form is simply a 
convection of the initial condition whereas the splitting method and the particle 
simulations agree on a different and more complex solution. 

These results show first that the MV model well describes the microscopic model 
in the dense regime. Secondly, that the correct formulation of the MV model is given 
by the limit of a conservative equation with a stiff relaxation term. 

The theoretical and numerical studies of the MV model highlight the specificity 
of non-conservative hyperbolic models with geometric constraints. More theoretical 
work is necessary in order to understand why the splitting method matches the 
microscopic model whereas the other methods do not. In particular, an extension 
of the theory developed in [8] to non-conservative relaxed models would be highly 
desirable. 

The outline of the paper is as follows: first, we present the Vicsek and MV models 
in section [2j Then, we analyze the MV model and give two different formulations 
of the model in section [31 We develop different numerical schemes based on these 
formulations and we use them to numerically solve different Riemann problems in 
section HJ Finally, we compare simulations of the microscopic model with those 
of the macroscopic system in the same situations in section Finally, we draw a 
conclusion. 
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2 Presentation of the Vicsek and Macroscopic 
Vicsek models 



At the particle level, the Vicsek model describes the motion of particles which tend 
to align with their neighbors. We denote by x k the position vector of the k th particle 
and by u k its velocity with a constant speed (|c<j fc | = 1). To simplify, we suppose 
that the particles move in a plane. Therefore x k G M 2 and u>k G S l . The Vicsek 
model at the microscopic level is given by the following equations (in dimensionless 
variables) : 

(2.D 

dujt = (Id — ujt ® uJk)(uJk dt + V2ddB t ), (2.2) 

where Id is the identity matrix and the symbol <E> denotes the tensor product of 
vectors, d is the intensity of noise, B t is the Brownian motion and is the direction 
of mean velocity around the k th particle defined by: 

&k = Tj-7 , Jk = Wj, (2.3) 

j,\xj-x k \<R 

where R defines the radius of the interaction region. Equation (12.21) expresses 
the tendency of particles to move in the same direction as their neighbors. The 
operator (Id — Uk <S> LOk) is the orthogonal projector onto the plane perpendicular to 
Uk- It ensures that the speed of particles remains constant. This model is already 
a modification of the original Vicsek model [29], which is a time-discrete algorithm. 

The Macroscopic Vicsek model (MV) describes the evolution of two macroscopic 
quantities: the density of particles p and the direction of the flow Q. The evolution 
of p and Q is governed by the following equations: 

d t p + V x -(c lP n) = 0, (2.4) 
p (d t tt + c 2 (tt ■ V x )fi) + A (Id - n ® tt)V x p = 0, (2.5) 
\tt\ = 1, (2.6) 

where ci, c 2 and A are some constants depending on the noise parameter d. The 
expressions of c\, c 2 and A are given in appendix |A] By contrast with the standard 
Euler system, the two convection coefficients C\ and c 2 are different. The other 
specificity of this model is the constraint = 1. The operator (Id — Q(g)f2) ensures 
that this constraint is propagated provided that it is true at the initial time. The 
passage from (123ft - (j2~3]) to (l2T4 >(j2^ - (pi)j) is detailed in [13]. We note that vortex 
configurations are special stationary solutions of this model in two dimensions (see 
appendix [B]). Up to our knowledge, this is the first swarming model which have 
such analytical solutions. 
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3 The Macroscopic Vicsek model 

3.1 Theoretical analysis of the macroscopic model 

To study model (l2.4D - (l2.5|) - (l2.6p . we first use the rescaling x' = xjc\. Then equa- 
tions (J23D-([23D-(J21D are wri tten: 



dtp + V x * ■ (pU) = 0, 

p (d t n + c'(n ■ Va./)fi) + a' (id - n ® n)v x <p = o, 

M = 1, 



(3.1) 
(3.2) 
(3.3) 



with d = oijc\ and A' = X/c\. In the sequel, we drop the primes for clarity. We 
refer to the appendix |A] for the computation of c and A and we just mention that 
we have (see figure fT71) : 



< c < 1 and A > 0, for all d > 0. 



(3.4) 



In two dimensions, we can use a parameterization of Q in polar coordinates: 
Q = (cos9,sm9) T . Therefore, equations f l3.ip - fl3.2l) can be rewritten as: 



dtp + d x (p cos 9) + d y (p sin 9) 
d t 9 + c cos 9d x 9 + c sin 9d y 9 + A 



0, 

sin 9 

P 



OxP H ~°vP 



(3.5) 
0. (3.6) 



In this section, we suppose that p and 9 are independent of y meaning that we are 
looking at waves which propagate in the x-direction. Under this assumption, the 
system reads: 



d t ( P e ^+A{p,9)d x ( P e ^ =0, 



with 



cos 9 —psin9 

_Asine ccog 



The characteristic velocities of this system are given by 



71,2 



c+ l)cos#± J(c- l) 2 cos 2 ^ + 4Asin 2 ^ 



(3.7) 
(3.8) 

(3.9) 



with 71 < 72. Therefore, the system is strictly hyperbolic. A possible choice of right 
eigenvectors is 



psin 9 
cos # — 71 



r 2 



c cos 6* — 72 

A sin 



(3.10) 



The two fields are genuinely nonlinear except at 9 
values of 7 p which satisfy: 



tan 2 6» 



J_ f ((c-l) 2 -4A)- 
4A [ (c+1) 2 



0, 9 = 7r and at the extrema 



(c 




Figure 1: The two eigenvalues 71 and 72 depending on 9 (d = 1 in this graph). 
For each curve, there exists a unique extremum (9± and 62) which corresponds to a 
degeneracy of the system. 

The Riemann invariants of the system (13. 7p are given by: 

I sin 5 

h[p ,e) = i°g/>-| — — V) ds, (3.11) 

h(p,e) = log,- f ^jM is. (3.12) 
Je A sm s 

The integral curve W\ and W2 starting from (pi, 9{) are given by: 

(f sin s A 
/ rr ds > ( 3 - 13 ) 
Je coss- 7l (s) y 

/ ccos s — 79 (s) , \ /„ „ ,s 

p 2 (0) = Piexp / . . m ; da . 3.14 

These are the rarefaction curves. To select the physically admissible rarefaction 
curve, we remark that 7 p must grow from the left to right states. The proofs of 
these elementary facts are omitted. The quantities 712 as functions of 9 are depicted 
in figure [TJ 
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3.2 A conservative form of the MV model in dimension 1 



For non-conservative systems, shock waves are not uniquely defined [21,23]. How- 
ever, in the present conservative formulation of the system can be found in 
dimension 1. Indeed, it is an easy matter to see that, if sin 9 ^ 0, system ( 13. 7p can 
be rewritten in conservative form: 



dt 



with: 



P 

W) 
Mo) 



d x 



p cos 9 
cf 2 (9) - Xlog(p) 



0. 



log 



tan ■ 



log 



sin^ 



cos 9 + 1 



log | sin 9 1 



(3.15) 

(3.16) 
(3.17) 



However, the functions f\ and f 2 are singular when sin 9 = which means that the 
conservative form is only valid as long as 9 stays away from 9 = 0. 

The conservative form (13.151) leads to the following Rankine-Hugoniot conditions 
for shock waves: two states (pi,9i) and (p r ,9 r ) are connected by a shock wave 
traveling at a constant speed s if they satisfy: 



Pr ~ Pi 
fl(9 r )-fl(9l) 



cf: 



2Wr 



p r COS 9 r — pi COS 9\ 

- c/ 2 - A log p r + A log pi 



(3.18) 



We can combine the two equations of the system (13.181) to eliminate the constant 
s and this leads to the following expression of the shock curve: 



{Pr ~ Pl)(cf 2 (9 r ) -Cf 2 



i) - Alogp r + A log pi) 

(p r COS# r - pi COs9i)(f 1 (9 r ) ~ /l (#/))■ 



(3.19) 



This equation must be numerically solved. The entropic part of the shock curve 
is determined by the requirement that 7 P must satisfy the Lax entropy condition. 
In figure [21 we give an example of a solution of a Riemann problem obtained by 
computing the intersection of the shock and rarefaction curves. 



3.3 The MV model as the relaxation limit of a conservative 
system 

We are going to prove that the MV model (I3.1l) - (l3.2p - (l3.3l) can be seen as the 
relaxation limit of a conservative hyperbolic model with a relaxation term. This 
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0.5 1 1.5 2 2.5 

P 



Figure 2: A solution of the Riemann problem with left and right states Ui and U r 
(solid line for shock waves and dotted line for rarefaction waves). In this example, 
the solution is given by two shock waves. 

link will be used later to build a new numerical scheme. More precisely, we introduce 
the relaxation model: 



In this model, the constraint \Q\ = 1 is replaced by a relaxation operator. Formally, 
in the limit e — > 0, we recover the constraint \Q\ = 1. 

Proposition 3.1 The relaxation model A3.2(jfy - fi3.2U\) converges to the MV model 
(I3.1j) - (l3..2j) -( l3.3)) as e goes to zero. 

Proof (formal). We define R £ = p £ (l — \Q £ \ 2 )Q £ . Suppose that as e goes to zero: 
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\p £ + v £ • (p £ n £ ) = o, 



(3.20) 



d t (p £ n £ ) + cv x ■ (p £ n £ ® n £ ) + x v xP £ 



- \n £ \ 2 )n £ . 



(3.21) 





(3.22) 



Then R £ 0, which generically implies that |f2°| 2 = 1 (except where p°f2° = 1 
which one assumes to be a negligible set). Therefore, we have: 



d t n° ■ n 







-Q° 



0. 



(3.23) 



Then since R £ x fi £ = 0, we have: 



{d t (p £ n £ ) + cV x ■ (p £ n £ ® tt £ ) + A V^p") x tt £ = 0. 



s 



and consequently when e — > 0: 

d t (p°Q°) + cV x ■ (p°Q° ® + A V x .p° = (3.24) 

for a real number a to be determined. Taking the scalar product of (13.241) with Q° 
and using (13.231) . we find: 

a = d t p° + cV x ■ {p°Q°) + \V x p° ■ n°. 
Using the conservation of mass (d t p° = — V x . • (p°f2 )), we finally have: 

a = (c - 1) V* • {p°Q°) + XV x p° ■ n°. 
Therefore, the relaxation term satisfies: 

-^R e = [( C - i)v x • (p°n°) + xv xP ° ■ + o(e). 

Inserting in (I3.20I) - (I3.21I) and taking the limit e — > 0, we recover the MV model 
(l3TT ]l - (p^ at the first order in e. 

□ 

Remark. As for the MV model, we can also analyze the hyperbolicity of the left 
hand side of fj3.20ri - fl3.2in . The eigenvalues are given by: 

7i = cu - \/A , 72 = cm , 73 = cu + VA, 

where u denotes the x-coordinate of Q and A = A — (c — c 2 )u 2 . The system is 



hyperbolic if and only if \u\ < y^2- As we can see in figure El for u 2 = 1, A is 
positive for any values of the noise parameter d. In particular, this implies that the 
relaxation model is hyperbolic for every \u\ < 1. 

4 Numerical simulations of the MV model 

4.1 Numerical schemes 

We propose four different numerical schemes to solve the MV model. The first two 
schemes originate from the discussions of the previous section, the two other one 
are based on the non-conservative form of the MV model. 

We use the following notations: we fix a uniform stencil (xj)j (with — Xi\ = Ax) 
and a time step At. We denote by £7™ = (p", #") the value of the mass and flux 
direction at the position Xj and at time nAt. 
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Figure 3: The quantity a/A/(c — c 2 ) depending on d. The relaxation model f)3.20p - 
f)3.20p is hyperbolic when the speed \Q\ is below this curve. At the limit e — > 0, 
\Q £ \ — > 1 and therefore the relaxation model is hyperbolic for any d in this limit. 



4.1.1 The conservative scheme 

Here we use the conservative form of the MV model ( 13. 151) : 

d t V + d x F(V) = 0, 



(4.25) 



with^ = (p,/i(^)) T andF(\/) = (pcos#, cf 2 (9)-X log(p)) T . We use a Roe method 
to discretize this equation: 



n+l 



vr Fi_ 
+ — 



At Ax 
where the intermediate flux Fi + is given by: 



~ FQQ + F{V l+1 ) 
F i+ - 

and A is the Jacobian of the flux F: 

A(V) = DF(V) = 



cos 6 —p sin 2 6 

-- ccos^ 
p 



(4.26) 



(4.27) 



(4.28) 



calculated at the mean value Vj 



Vi+V i+1 



As mentioned earlier, the conservative form is only valid when 9 does not cross a 
singularity 9 = or 6 = tt (i.e. sin 9 = 0). Nevertheless, numerically we can still use 
the formulation (I4.26P when 9 changes sign. Moreover, since f\ is an even function, 
this only gives |# n+1 |. To determine the sign of 9 n+1 , we use an auxiliary value 9 
which we update with the upwind scheme (14.321) . The sign of 9 is then determined 
using the sign of 9. 
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4.1.2 The splitting method 

The next scheme uses the relaxation model (I3.20p - fl3.2ip . The idea is to split the 
relaxation model in two parts, first the conservative part: 

d tP + V, • (pQ) = 0, 

d t (pfl) + cV x ■ (pQ ® Q) + A V x p = 0. 1 } 

and then the relaxation part: 

dtP =0, 

This system reduces to: dtVL = -(1 — \Q\ 2 )Q. Since this equation only changes the 
vector field Q in norm (i.e. d t fl ■ fi -1 = 0), we can once again reduce this equation 
to: 

^d t \n\ 2 = l(i - \n\ 2 )\n\ 2 . (4.3i) 

Equation flOTD can be explicitly solved: \VL\ 2 = (l+C e~ 2 / £t y l withC = - 1 

We indeed take the limit e — > of this expression and replace the relation term by 
a mere normalization: f2 — > Q/\Q\ . 

The conservative part is solved by a Roe method with a Roe matrix computed 
following [23] page 156. 

4.1.3 Non-conservative schemes 

We present two other numerical schemes based on the non-conservative formulation 
of the MV model. 

(i) Upwind scheme 

The method consists to update the value of [/" with the formula: 

rm+l _ Tin /Tjn_Tjn \ /jjn _Tjn\ 

"i . A+ ( Ui u *~i \ +A -( u * + i u _a ] =0; (432) 



At V At / \ At 

where A + and A~ are (respectively) the positive and negative part of A, defined 
such that A = A + — A~ and \A\ = A + + A~ and A + , A~ are computed using an 
explicit diagonalization of A. 
(ii) Semi-conservative scheme 

One of the problem with the upwind scheme is that it does not conserve the 
total mass (j x p(x) dx). In order to keep this quantity constant in time, we use the 
equation of conservation of mass A3. If) in a conservative form: 

d t p + d x H(p : 6)=0, (4.33) 
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with H(p,9) = pcos9. Therefore, a conservative numerical scheme associated with 
this equation would be: 



p F\-ti + H W- H "I* = o, (4.34) 
At Ax 

where H i+ is the numerical estimation of the flux H at the interface between Xi 
and Xi + i. To estimate numerically this flux, we use the following formula with 

Ui = {Pi,9i): 

^ fTjn _ Tjn \ 

H t+1/2 = H(U l+l/2 ) - \A\ p ( m 2 - j , (4.35) 

where the intermediate value is given by fj+1/2 = ' 2 1+1 and \A\ is the first line 
of the absolute value of A. 

For the estimation of the angle 9, we keep the same scheme as for the upwind 
scheme. This numerical scheme uses one conservative equation (for the mass p) 
and a non-conservative equation (for the angle 9). It is thus referred to as the 
semi-conservative scheme. 



4.2 Numerical simulations 

To compare the various numerical schemes, we use a Riemann problem as initial 
condition. We choose solutions which consist of a rarefaction wave (figure H]) or a 
single shock wave (figures [MS]) . 

We take the following parameters: d=l, the length of the domain is 10 units and 
the discontinuity for the Riemann problem is at x = 5 (the middle of the domain). 
The simulations are run during two time units with a time step At = 2.10 -2 and a 
space step Ax = 5.10 -2 . For these values, the Courant number (C n ) is 0.778. We 
use homogeneous Neumann conditions as boundary conditions. 

For the rarefaction wave, we take: 

{Ph 0i) = (2, 1.7) , (p n 9 r ) = (1.12, 0.60). (4.36) 

All the numerical schemes capture well the theoretical solution (see figure H]). 
For the shock wave, we choose: 

{p h B{) = (1, 1.05) , (p n 9 r ) = (1.432 ,1.7). (4.37) 

For these values, the shock speed computed with (13.191) is s = —1.585. The re- 
sults of the numerical simulations using the four schemes are given in figure The 
numerical solutions are in accordance with the theoretical solution given by the con- 
servative formulation for all the numerical schemes. Nevertheless, the conservative 
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scheme is in better accordance with this solution. For the other schemes, the shock 
speed differs slightly. 

A second example of a shock wave is computed using the following initial condition: 

{p h 0,) = (1,0.314) , (p PJ r ) = (2,1.54). (4.38) 

The solutions given by the 4 numerical schemes are very different. Only the con- 
servative method is in agreement with the solution given by the conservative for- 
mulation. But the conservative formulation is not necessary the right one. Indeed, 
in the next section, particle simulations show that the right solution is not given by 
the conservative formulation but rather by the splitting method. 




Figure 4: The theoretical solution of the Riemann problem (I4.36P given by a rar- 
efaction curve (solid line) and the numerical solutions (points), p (blue) and cos# 
(green) as functions of space. The simulations are run during 2 time units, with a 
time step At = 2.1(T 2 and a space step Ax = 5.10" 2 (CFL=.778). 
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ure 5: Theoretical and numerical solutions of the Riemann problem ( 14. 3 7ft 
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ure 6: Theoretical and numerical solutions of the Riemann problem (I4.38I) 
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5 The microscopic versus macroscopic Vicsek 
models 



5.1 Local equilibrium 



In this part, we would like to validate the macroscopic Vicsek model by the simu- 
lation of the microscopic Vicsek. The macroscopic model relies on the fact that the 
particle distribution function is at local equilibrium given by a Von Mises distribu- 
tion M u (see [13]): 

M n (u) = Cexp(^J (5.1) 

where C is set by the normalization condition^. The goal of this section is to show 
numerically that the particle distribution of the microscopic Vicsek model is close 
in certain regimes to this Von Mises distribution. 

To this aim, in appendix Owe propose a numerical scheme to solve system ( 12.21) . 
The setting for our particle simulations is as follows: we consider a square box with 
periodic boundary conditions. As initial condition for the position Xi, we choose a 
uniform random distribution in space. The velocity is initially distributed according 
to a uniform distribution on the unit circle. 

During the simulation, we compute the empirical distribution of the velocity 
direction 9 and of the mean velocity Q of particles. We then compare this empirical 
distribution with its theoretical distribution Mq(6) given by (15.11) . In figured we 
give an example of a comparison between the distribution of the velocity direction 
9 and the theoretical distribution Mq predicted by the theory. 

Since the distribution of velocity converges, we have a theoretical value for the 
mean velocity. We denote by tp^ the mean velocity of particles and ip the theoretical 
value given by the stationary distribution: 



1 

N 



N 
k=l 



LO Ma(uj) duj 



(5.2) 



At least locally in x, we have that (p^ <f- In figure we compare the two distri- 
butions for different values of the noise d and we can see that the two distributions 
are in good agreement. We also observe a smooth transition from order (<p « 1) to 
disorder (</? << 1) as it has been measured in the original Vicsek model [29]. 

The situation is different when we look at a larger system. We still have conver- 
gence of the velocity distribution of particles to a local equilibrium p(x) Mj 



n(x) 



but the mean direction £l(x) now depends on x. Therefore the mean velocity of the 



Explicitly given by C 1 = 2ir Io(d x ) where Iq is the modified Bcssel function of order 
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Figure 7: Left figure: the distribution of velocity direction 9 (with uj = (cos 9, sin#)) 
compared with its theoretical distribution after 6 time units of simulation. Right 
figure: the corresponding particle simulation. Parameters of the simulation: Lx = 
1, Ly = 1 (domain size), number of particles N = 500, £ — 1/4, R = .5, d = .2, 
At = 2.1Q- 3 . 




Figure 8: The mean velocity tp (15. 2p for different values of d. Parameters of the 
simulation: Lx = 1, Ly = 1 (domain size), number of particles iV = 200, radius 
of interaction R = .5, At = .02 unit time, the simulations are run during 180 unit 
time. 
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particles in all the domain differs from the expected theoretical value (I5.ip - (I5.2I) . 
We illustrate this phenomena in figure EJ we fix the density of particles and we 
increase the size of the box. As we can observe, the mean velocity <fiN (15. 2p has a 
smaller value when the size of the box increases. This phenomena has been pre- 
viously observed in [7]. The mean velocity </?jv can also differ from the expected 
theoretical value cp (I5.ip - (I5.2I) when the density of particles is low. In figure [TUl we 
fix the size of the box (L = 10) and we increase the density of particles (the density 
is given by the number of particles inside the circle of interaction). At low density, 
the mean velocity tp^ is much more smaller than the theoretical prediction <p. But 
as the density of particles increases, the mean velocity ip^ grows (see also [29]) and 
moreover ip^ converges to <p. Because of that, a dense regime of particles has to be 
used in the following in order to numerically compare the microscopic model with 
the MV model. 



5.2 Microscopic versus Macroscopic dynamics 

We now compare the evolution of the two macroscopic quantities p and Q for the two 
models. We have seen that the different schemes applied to the macroscopic equa- 
tion can give different solutions (see figure [6]). Therefore, we expect that particle 
simulations will indicate what is the physically relevant solution of the macroscopic 
equation. 

We first briefly explain how we proceed to run the particle simulations of a Rie- 
mann problem (see also appendix O) . First, we have to choose a left state (pi,9i), 
a right state (p r , 9 r ) and the noise parameter d. Then we distribute a proportion 
^ of particles uniformly in the interval [0, 5] and the remaining particles uni- 
formly in the interval [5,10]. Then, we generate velocity distribution 9 for the 
particles according to the distribution Mq (15. ip with Q t = (cos0/, sin^) T on the 
left side and Q r = (cos0 r , sin0 r ) T on the right side. We use the numerical scheme 
given in appendix O to generate particle trajectories. To make the computation 
simpler, we choose periodic boundary conditions. Therefore the number of parti- 
cles is conserved. As a consequence, there are two Riemann problems corresponding 
to discontinuities at x = 5 and at x = or 10 (which is the same by periodicity). 
We use a particle-in-cell method [15,19] to estimate the two macroscopic quantities: 
the density p and the direction of the flux Q (which gives 9) . In order to reduce the 
noise due to the finite number of particles, we take a mean over several simulations 
to estimate the density p and 9 (10 simulations in our examples). 

In figure [TTJ we show a numerical solution for the following Riemann problem: 

{ Ph 0,) = (1, 1.5) , (p r , 9 r ) = (2, 1.83) , d = 0.2 (5.3) 

using particle simulations and the macroscopic equation. We represent the solutions 
in a 2D representation. Since the initial condition is such that the density p and 
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Figure 9: The mean velocity tp (15. 2p for different values of d. We use different do- 
main sizes and we keep the same density of particles. As the domain size increases, 
the total flux tp decreases which means that particles are less aligned globally. Pa- 
rameters of the simulations: L = 1, 2, 5, 10 (domain size), number of particles 
N = 200, 800, 5000, 20000, radius of interaction R = .5, At = .02 time units, the 
simulations are run during 180 time units. 
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Figure 10: The mean velocity tp (15. 2p for different values of d. We change the 
density of particles (given by the mean number of neighbors in unit of radius of 
interaction). When we increase the mean number of neighbors, particles are more 
aligned. Parameters of the simulations: L = 10 (domain size), number of particles 
N = 254, 1273, 6366, 12732, and 25464, radius of interaction R = .5, At = .02 time 
units, the simulations are run during 180 time units. 
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the direction 9 are independent of the y-direction, we only represent p and 9 along 
the x-axis in the following figures. 

In figure [T2J, we represent the two solutions (the particle and the macroscopic 
one) with only a dependence in the x-direction. Three quantities are represented: 
the density (blue), the flux direction 9 (green) and the variance of the angle dis- 
tribution (red). The macroscopic model supposes that the variance of 9 should be 
constant everywhere. Nevertheless, we can see that the variance is larger in regions 
where the density is lower. For p and 9, we see clearly the propagation of a shock in 
the middle of the domain and a rarefaction at the boundary. The CPU time for one 
numerical solution at the particle level is about 140 seconds with the parameters 
given in figure [12j For the macroscopic equation, the CPU time is about 0.1 second 
which represents a cost reduction of three orders of magnitude compared with the 
particle simulations. Since we have to take a mean over many particle simulations, 
the cost reduction is even larger. 

In figure [131 we use t ne same Riemann problem to set the initial position as 
in figure [6] both with d = 1 (14.381) . We use a larger domain in x (L = 20 space 
units) in order to avoid the effect of the periodic boundary condition. The upwind 
scheme and the semi-conservative method are clearly not in accordance with the 
particle simulations. Moreover, the splitting method is in better agreement with 
the particle simulation since the shock speed is closer to the values given by the 
particle simulations than that predicted by the conservative scheme. 

Finally, our last simulations concern a contact discontinuity. We simply initialize 
with: 

(fl,0,) = (l,l) , (p rj r ) = (l,-l) , d = 0.2, (5.4) 

i.e. we reflect the angle with respect to the x-axis across the middle point x = 5. A 
natural solution for this problem is the contact discontinuity propagating at speed 
ccos(l): 

p(t,x) = l , 9(t,x) = 9 (x - ccos(l)t), (5.5) 

with 9 (x) = — 1 when x < 5 and 9 (x) = 1 when x > 5. This is the solution pro- 
vided by the conservative scheme (figure [T41) . But surprisingly, the splitting method 
and the particle simulation agree on a different solution. Indeed, the solutions given 
by the particles and the splitting method are in fairly good agreement with each 
other, which seems to indicate that the "physical solution" to the contact problem 
(15 .4p is not given by the conservative formulation (I5.5P but by a much more complex 
profile. The constraint of unit speed drastically changes the profile of the solution 
compared with what would be found for a standard system of conservative laws. 
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Density (rho) 



Density (rho) 




Figure 11: The particle density in space p computed with particle simulations (left) 
and the macroscopic equations (right). We initialize with a Riemann problem (15.31) . 
Numerical parameters for the particle simulations: N = 2.10 6 particles, At = .01, 
e = 1/10, R = .5, Lx = Ly = 10, we take a mean over 10 computations. Numerical 
parameters for the macroscopic model: At = .01, Ax = .025 (CFL=0.416), we use 
the splitting method. The simulations are run during 2 time units. 




Figure 12: The solution of the Riemann problem (15.31) with d = .2 computed with 
the splitting method (solid line) and with particle simulations (dots). In blue, we 
represent the density p, in green the flux direction 9 and in red the variance of the 
velocity direction. The parameters are the same as in figure [TTJ We only change 
the representation of the solution (ID-representation). 
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Figure 13: The solutions of the Riemann problem (I4.38P with d — 1 computed 
using the macroscopic model and particle simulations of the microscopic model (see 
figure [T2|) . Numerical parameters for the macroscopic model: At = .01, Ax = 
.025 (CFL=0.778). Numerical parameters for the particle simulation: N = 2.10 6 
particles, At = .02, e = .1, R = .5, Lx = 20 and Ly = 1. We take a mean over 50 
simulations. The simulations are run during 6 time units. Since d — 1, fluctuations 
are higher (see figure [9]), we have to increase the density of particles to reduce this 
effect. 
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Figure 14: The solution of the Riemann problem f ]5.4f) computed with the conser- 
vative method (top), the splitting method (down) and with particle simulations 
(dots). Numerical parameters for the macroscopic model: At = .01, Ax = .025 
(CFL=0.416). Numerical parameters for the particle simulations: N = 10 6 parti- 
cles, At = .01, e = 1/10, R = .5, Lx = 10, Ly = 1. We take a mean over 100 
simulations. The simulations are run during 2 time units. 
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6 Conclusion 



In this work, we have numerically studied both the microscopic Vicsek model and 
its macroscopic version [13]. Due to the geometric constraint that the velocity 
should be of norm one, the standard theory of hyperbolic systems is not applicable. 
Therefore, we have proposed several numerical schemes to solve it. By comparing 
the numerical simulations of the microscopic and macroscopic equations, it appears 
that the scheme based on a relaxation formulation of the macroscopic model, used in 
conjunction with a splitting method is in good agreement with particle simulations. 
The other schemes do not show a similar good agreement. In particular, with an 
initial condition given by a contact discontinuity, the splitting method and the 
microscopic model provide a similar solution which turn to be much more complex 
than what we could be expected. 

These results confirm the relevance of the macroscopic Vicsek model. Since the 
CPU time is much lower with the macroscopic equation, the macroscopic Vicsek 
model is an effective tool to simulate the Vicsek dynamics in a dense regime of 
particles. 

Many problems are still open concerning the macroscopic Vicsek model. We 
have seen that the splitting method gives results which are in accordance with par- 
ticle simulations. But, we have to understand why this particular scheme captures 
well the particle dynamics better than the other schemes. Since the macroscopic 
equation has original characteristics, this question is challenging. Another point 
concerns the particle simulations. We have seen that the particles density has a 
strong effect on the variance of the velocity distribution. When the density is low, 
the variance is larger. The macroscopic equation does not capture this effect since 
the variance of the distribution is constant. Works in progress aims at taking into 
account this density effect. 
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A The coefficients ci, C2 and A 

The analytical expression of the coefficient c\ involved the distribution of the local 
equilibrium Mq (15. ip . The two other coefficients (c 2 and A) involve also the solution 
g of the following elliptic equation [13]: 

- (1 - x 2 ) d x [e x/d (l - x 2 )d x g] + e x/d g = -(1 - x 2 ) 3/ V /d , (A.l) 

on the interval x G (—1, 1). 

If we define the function h = ^f=s an d M(x) = e3, these macroscopic coeffi- 
cients can be written as: 

f n 7r cos0M(cos0)sin0d0 , s 

Cl = <cos0> |m= Jo r : ... , (A.2) 

J M (cos 0) sin d0 

cos sin 2 /i(cos 0)M(cos 0) sin 6 d6 

c 2 = < cos0 > \ sin 20 hM = r , 2 — — — , — , (A. 3) 

J sm 0/i(cos0)M (cos 0)sm 0c(0 

X — d. (A.4) 

In the above expressions, we can see that < c\, c% < 1. 

Now we are going to explore two asymptotics of g when the parameter d is small 
or large. 

Lemma A.l Let g be the solution of equation (fAJJ) . We iave tie asymptotics: 



d->0 , 

~ d 



asjn(x i ) — — 



o(d), (A.5) 



Proof (formal). Introducing the Hilbert space: 

V = {g | (1 - ^ 2 )- l/2 g e L 2 (-l, 1), (1 - ^ 2 d,g e L 2 (-l, 1)} 

we have already seen in [13] that there exists a unique solution g of (lA.ip . Moreover 
this solution is negative. 

To derive the asymptotic behavior of g depending on d, we first develop (1A. if) : 

d x [(l - x 2 )d x g] + (1 - x 2 ) - d 8 x g - = (1 - x 2 ) 1 / 2 . (A.7) 

When d — ► 0, we have: 

d x g = 

on the interval [—1 + e, 1 — e] for all e > 0. Since g belongs to V, we also have the 
boundary condition g( — l) = g(l) = 0, so g converges to when d — > 0. 
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To derive the next order of convergence in the limit d — > 0, we normalize g with 
g = dg, which gives: 

d d x [(l - x 2 )d x g] + (1 - x 2 )d x g - d -^— 2 = (1 - x 2 ) 1 / 2 . (A.8) 

In the limit d — > 0, we deduce that: 

{I - x 2 )d x g = {I - x 2 ) l '\ (A.9) 

which has an explicit solution: g = asin(x) + c. Since g < 0, we have c < — |. 
Numerically, we find that c = — | but the proof is still open. This formally proves 

When d — > +oo, ( 1A.7I) gives: 

d x [(l - x 2 )d x g ] - = (1 - x 2 ) 1 / 2 . (A.10) 

1 — or 



A simple calculation shows that go = — | yl — x 2 is a solution of flA.lOj) . 

To derive the next order of convergence, we look at the difference v = d(g — g ), 

which satisfies (see f ]A.7j) and (lA.lOj) ): 

d x [(l - x 2 )d x v] + (1 - x 2 ) - d d x v - = -(1 - x 2 ) d x g . 

In the limit d — > +oo, t> satisfies: 



d x [(l - x 2 )d x v] - j = --xvT^. (A.ll) 

1 — x z 2 

A simple calculation shows that v = y^xVl — x 2 is solution of ( 1A.11I) . Therefore 
we formally have the expression (1A.6I) in the proposition. 

□ 

In figure [TS1 we compute numerically the function g (lA.lj) . We use a finite 
element method with a space step Ax = 10 -3 . The two asymptotics of g when 
d —* and d — » +oo are computed in figure | 



Proposition A. 2 The two coefficients C\ and C2 defined (resp.) by the equations 
(\A.2\) and hA.3\) satisfy the following asymptotics: 



ci ~ l-d + 0(d 2 ), (A.12) 
3d V d 2 



<* ~ i + O^ , (A.13) 



C2 6^ + °U)- 
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Proof. We have an explicit expression for the coefficient c\ using the change of 
unknowns x = cos(0): 

ci = coth - d, (A. 15) 

where coth(s) = e ^ e _^ . The expressions of ( 1A.12[) and (1A.13I) are simply deduced 
by a Taylor expansion of the last expression. 

For the coefficient C2, we insert the development of g (1A.6j) in expression (1A.3I) . 

□ 

Remark. The behavior of c 2 when d — > is more difficult to analyze. The 
density probability sin 6hM used in formula (1A.3I) becomes singular in this limit. 
Nevertheless, due to the expression of Mq, the density converges to a Dirac delta at 
which explains why C2 d ~° 1. To capture the next order of convergence, we need 
to find the second order correction of g in the limit d — > which is not available. 
However, numerically we find that: 

c 2 1 -2d + o(d). 

In figure [T71 we numerically compute the coefficients c^/cx, \/c\ and their asymp- 
totics. 




-1 -0.5 0.5 1 

x 



Figure 15: The numerical solution g of (1A.1I) for different values of the parameter d. 
We have the following asymptotics (see lemma |A~TT) : g ^> and g — — x 2 . 
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Figure 16: Left figure: the first correction g\ of g when d — > 0. The red curve is 
the theoretical asymptotic limit: g± — g/d asin(x) — | (see lemma lATTT) . Right 
figure: the first correction gi of g when d — > +oo. The red curve is the theoretical 

asymptotic: g\ = d(g + \ y/1 — x 2 ) y^Vl — x 2 (see lemma lATTT) . 
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Figure 17: The ratio c<ijc\ and \ jc\ (solid lines) and their two asymptotics (dashed 
lines) (see proposition IA. 21) computed with Ax = 10 -3 . 



27 



B Special solution of the MV model 

In this appendix, a vortex configuration is exhibited as a stationary solution of the 
MV model fl2T4l) - fl2~5|) - fl2T6l) in dimension 2. A stationary state of the MV model 
has to satisfy: 

v,-(p^) = o, 

c(n- v x )n + \(u-n®n)?f = o. 1 ' 

Introducing the polar coordinates, p(r, 9), Q(r, 9) = f r (r, 9)e r +fg(r, 9) eg, where e r = 
(cos#,sin#) T and eg = (— sin 9, cos 9) T , we are able to formulate the proposition: 

Proposition B.l The following initial condition is a stationary state of the MV 
model ( fB.Ij) : 

p(r) = Cr c / A , n = e e , (B.2) 

where C is a constant. 



Proof. With the expression of p and f2 given by ( 1B.2j) . the divergence of the mass 
is zero and the gradient of p is orthogonal to fl, therefore the system (IB.lj) reduces 
to: 

c(fi • V x )fi + = 0, (B.3) 

or in polar coordinates: 

1 p'(r) 
c-dgeg + \^H-e r = 0. 
r p(r) 

Since 9eee = — e r , we can easily check that the solution of this equation is given by 
p(r)=Cr c / x . □ 



C Numerical schemes for particle simulations 

In the limit e — > 0, an explicit Euler method for the differential system fl2.1|) - fl2.2p 
imposes a restriction time step condition of ^At < 1. Therefore, we develop an 
implicit scheme for this system. The idea is to go back to the original Vicsek model 
(see [13]). We use the formulation: 

, ,n+l .n 

Af = (M - u n+1/2 <8> u n+1/2 ) (u n - u n ) (C. 1) 

where uj n+l ^ 2 = r^j^+ri and u n is the average velocity (12.31) . When At = 1, we 
recover exactly the orig inal Vicsek model [29]. (lOH can in fact be solved explicitly. 
First, we have to remember that uo n+l belongs to the unit circle (i.e. |a> n+1 | = 1). 
Then we use that u n+1 — uj 11 is the orthogonal projection of (u n — uo n )At on the 
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orthogonal plan of uj n+l / 2 . Therefore u n+1 and u n are on the circle C with center 

(see figure [TSi) . This fully defines u n+1 



B = uj n + — — ^ ^ At and radius 



{Q n -LJ n )At 

2 



since u n and uj n+l are the two intersection points of the unit circle and the circle 
C. Denoting 9 the angle of the unit vector u, we easily check that we have in terms 
of angles: 

e n+l = r + 2 (^TB). 

To take into account the effect of the noise, we simply add a random variable: 



e n+i = e n + 2 (u n , B) + V2dAte n 
where e n is a standard normal distribution independent of 9 n . 



(C.2) 




{Q n - u n )At 



c ,-■ 



Figure 18: Illustration of the geometric method to solve explicitly equation fIC.lj) 



Algorithm used to solve a Riemann problem with particles. 

1. Choose a Riemann problem (pi, 8i) and (p r , 8 r ). 

2. Initiate TV particles (x^, ujk)k=i. .n according to the distributions piM^ and 
p r M Qr . 



3. Let evolve the particles in time using the time-discretization (10.2[) of equation 

4. Compute the mass p and the direction of the flux Q using Particle-In-Cell 
method [19] in order to compare the simulation with the one of the MV 
model. 
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